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FOREWORD 


The  work  reported  herein  was  sponsored  by  the  Arnold  Engineering 
Development  Center  (AEDC),  Air  Force  Systems  Command  (AFSC), 
under  Program  Element  62405214.,  Project  6951,  Task  695102. 

The  results  of  the  research  presented  were  obtained  by  ARO,  Inc. 

(a  subsidiary  of  Sverdrup  and  Parcel,  Inc. ),  contract  operator  of  the 
AEDC,  AFSC,  Arnold  Air  Force  Station,  Tennessee,  under  Contract 
AF40{600)-1200.  The  research  was  performed  from  January  until  June, 
1965,  under  ARO  Project  No.  RW2408,  and  the  manuscript  was  submitted 
for  publication  on  November  17,  1965. 

The  solution  of  the  mathematical  problem  to  which  this  report  is 
addressed  was  stimulated  by  a  requirement  to  solve  a  set  of  ordinary, 
first-order,  nonlinear,  stiff  differential  equations  which  describe 
chemical  reaction  rates  applicable  to  hydrogen-air  reactions.  The 
authors'.wish  to  thank  R.  P.  Rhodes,  Jr.,  who  proposed  the  problem 
which  led  to  the  development  of  this  technique,  and  who  worked  jointly 
with  the  authors  in  solving  the  proposed  problem. 

This  technical  report  has  been  reviewed  and  is  approved. 


Marion  L.  Laster 
Propulsion  Division 
DCS/Research 


Donald  R.  Eastman,  Jr. 
DCS/Research 
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ABSTRACT 


A  method  is  presented  for  numerically  integrating  a  system  of 
stiff,  first-order  differential  equations.  This  method  is  based  on 
transforming  the  set  of  dependent  variables  so  that  the  resulting  sys¬ 
tem  will  not  be  stiff;  the  transformed  system  is  then  integrated  by  the 
Eunge-Kutta  method.  The  resulting  procedure  is  often  appreciably 
faster  than  classical  methods  in  that  a  much  larger  step  size  is  allow¬ 
able  with  nominal  increase  in  step  computation  time.  Applications  and 
results  are  discussed  for  systems  of  various  order,  including  a  system 
of  six  chemical  rate  equations. 
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SECTION  I 
INTRODUCTION 


Although  systems  of  diffei  ential  equations  having  the  property  of 
being  stiff  have  been  used  as  mathematical  models  of  certain  physical 
phenomena  almost  since  the  invention  of  the  calculus,  only  in  recent 
years  have  such  systems  been  categorized  according  to  this  property 
and  been  formally  recognized  as  a  class  of  equations,  the  numerical 
solution  of  which  often  eludes  the  computer  when  sought  by  classical 
means.  Since  1951,  various  works  have  been  published  on  the  subject 
(Refs.  1  through  5). 

The  primary  reason  that  equations  in  this  class  present  such  dif¬ 
ficulty  is  that  an  exceedingly  small  integration  step  size  is  sometimes 
required,  making  classical  integration  techniques  impractical  even  on 
the  most  sophisticated  computing  machines. 

The  object  of  this  presentation  is  to  formulate  a  method  that  can 
be  used  efficiently  on  a  high  speed  digital  computer  to  obtain  the  solu¬ 
tion  of  a  system  of  first-order,  ordinary,  stiff  differential  equations. 
The  method  described  is  straightforward  and  practical  when  applied  to 
stiff  equations  resulting  from  many  physical  problems. 

The  above  objective  may  be  realized  by  transforming  the  dependent 
variables  in  such  a  manner  that  the  resulting  system  of  equations  will 
not  be  stiff  in  a  neighborhood  of  the  value  of  the  independent  variable  at 
which  the  transformation  was  made.  The  Runge-Kutta  method  will  be 
used  to  integrate  the  resulting  system. 


SECTION  II 

DEFINITIONS  AND  NOMENCLATURE 


When  presenting  a  similar  technique,  some  authors  develop  their 
theory  exclusively  for  a  single  equation  with  one  dependent  variable, 
and  later  state  that  the  extension  of  their  method  to  a  system  of  such 
equations  is obvious".  It  is  the  contention  of  the  authors  of  this 
presentation  that  any  such  extension  is  almost  always  |ambiguous.  It 
is  also  believed  that  with  little  or  no  sacrifice  to  clarity,  the  method 
can  be  developed  in  vector  notation,  leaving  no  doubt  concerning  how 
the  method  should  be  applied  to  the  general  case.  Following  these 
convictions,  the  nomenclature  introduced  below  Will  be  used  throughout 
this  document: 
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1.  Roman  characters  will  be  used  to  denote  real  scalar 
quantities. 

2.  Greek  characters  will  denote  numbers  which  are,  in 
general,  complex  scalars. 

—fc 

3.  An  (-.)  will  indicate  an  n  x  1  column  vector,  and  <f> 
denotes  the  null  vector. 

4.  An  (-)  will  denote  an  n  x  n  matrix;  ^  and  I  denote  the 
null  matrix  and  identity  matrix,  respectively. 

5.  The  notation  Re  (A)  will  mean  the  real  part  of  the 
complex  number,  A  . 

Now  consider  the  system  of  equations 

=  Fj  (x,  y,,  y3,  •  •  •  ,  yn)  ;  yi  (x0)  =  yi„  (1) 

for  i  =  1  n  .  More  conveniently,  Eq.  (1)  can  be  written 

=  F  (x,y)  ,  y  (\c)  =  y„  (2) 

where 


Yi 

X 

_ 1 

Yt 

F2  (x,  y  ) 

’ 

■+  ^ 

. 

F  (x  ,  y  )  = 

, 

yn 

Fn  (*,  y) 

In  order  to  get  Eq.  (2)  in  an  applicable  form,  define  the  matrix 

<3F, 


“  f  (x,y)  = 


’  <?Ft  , 

ar. 


OF, 


P-(x.y)  ^(x,y) 


&Yi 


-P-2  (x,  y )  (xj  y ) 

°vi  uYi 


(x,y) 


dF, 

w«(x-y] 


dFr 

<hn 


(x,y) 


(3) 
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It  will  be  assumed  throughout  that  7  (x,  y)  is  nonsingular,  this 
method  being  inapplicable  at  any  point  on  tine  solution  of  Eq.  (2}  for 
which  this  is  not  true.  Equation  (2)  can  now  be  written 

ii-  +  f(x,y)y=g(x,y)  ,  y  ( *0 )  =  y0 
dx 

where 

g  (x, y)  -  F (x, y)  +  f  (x,  y)  y 

The  form  of  Eq.  (4)  is  now  such  that  the  definition  of  a  stiff  differ¬ 
ential  equation  can  be  given.  Definition:  Let  x  =  x,  be  a  value  of  the 
independent  variable  in  the  region  of  interest,  and  let  A,,  •  ,  A„ 

be  the  eigenvalues  of  f  (xlty  (xj  .  If 

Re  (Aj)  >  >  0 

for  any  i,  then  Eq.  (4)  is  said  to  be  stiff  at  x  =  x,  . 


(4) 

(5) 


SECTION  III 

DERIVATION  OF  THE  TRANSFORMATION 


It  was  mentioned  previously  that  the  dependent  variable  in  Eq.  (4), 
y  (x0)  ,  will  be  transformed  so  that  if  |  x  -  x„  1  £  h  ,  then  the  transformed 
equation  will  not  be  stiff.  For  the  remainder  of  this  paper,  h-is  under¬ 
stood  to  be  the  integration  step  size  for  the  Hunge-Kutta  method  (fourth- 
order)  on  the  transformed  equations.  The  maximum  allowable  value  of 
h|will,  of  course,  depend  on  the  exact  nature  of  the  original  differential 
equation;  the  transformation  has  allowed  the  value  of  h  to  be  significantly 
increased  in  all  test  cases  tried  so  far. 


It  is  pointed  out  again  that  the  validity  of  the  method  depends  on  the 
nonsingularity  of  f  (x,  y). 


It  is  convenient  to  define  at  this  time  the  matrix 


u  (x)  =  I  -  f0  (x  “  x0)  *■ 


f0  (*  ~  *p)  "fg  (*  -  X0) 


0 '  2  ! 
which  is  sometimes  referred  to  as 


3  ! 


(6) 


u  (x)  =  exp  (-  f0  (x  -  x0)) 


for  obvious  reasons. 
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The  following  properties  of  the  infinite  matrix  series,  Eq.  (G),  can 
be  readily  verified  independent  of  fa  and  (x  -  xa) ; 

1.  The  series  converges 

2.  u  (x)  is  nonsingular 

3 .  u  (x)  !a  =  7a  5  (x) 

Note  also  that 

(?) 

(8) 

(9) 


d  +  fa  u  (x)  =  0  .  u  (x0)  =  f 

Now,  let  z  (x)  be  a  solution  of 


Hence 

z  =  5  (x)  {ya  -  f“l  g€)  +  g# 


Defining 

w  (x)  =  y  (x)  -  z  (x)  one  gets,  by  subtracting  Eq,  (8)  from  Eq.  (4) 


where 


~~dx  '  +  f  ( x *  y)  w  (x)  ~  V  (x,  y)  ;  w0  -  ^ 


v  (x,y)  =  g(x,y)  -  go  +  [r„  -  f(x,y)]  z{x) 


(10) 


(11) 


Premultiplying  Eq.  (10)  by  u  (x)  and  post  multiplying  Eq.  (7)  by  w(x) 


gives 


and 


u  (x)  +uf(x,y)w  =  «v 

dx 


w  +  lQ  3  (x)  w  = 


Subtracting, 


U  — w  +  u  r  ( x ,  y )  w  -  f0  u  (x)  w  =  uv 


(iTl)  (u  wj  +  (u  l)^ur(x,y)w-f0u(x)wJ  =  {u  ' )  u 

(u-1w)+u  ^  f"  ( x ,  y )  -  u  lfa  u  J 


—  —  I  -» 

W  =  II  V 


(12) 
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Defining 


-  —  i  -* 

t  =  u  w 


m  =  u  ~ 1  F  (x ,  y )  —  frJu 


Eq.  (12)  becomes 


d  t 
dx 


+  m  t 


=  «*"’  V  ;  t  (x0)  -  4> 


(13) 


Equation  (13)  is  the  desired  transformed  equation,  the  relation  be¬ 
tween  i  and  y  being 


y  =  u  (x)  t  (x*l  ■*-  ffl 


1  r  -  sw 
_ 


F0  +  y„ 


Setting 


G  (x, t  )  =  u  '  v  -  m  t  =  u  *~t  j^F  (x,  y)  -  F0  +  f„  (y  -  yo)J 
Eq.  (13)  can  be  written  as 


J-L 

dx 


«  G  ( x  ,  t  )  ;  t  ( Xc  )  *=  tj> 


An  interesting  consequence  of  the  transformation  is  observed  if 
Eq.  (2)  is  linear  with  constant  coefficients.  In  that  case 

f  ( x ,  v )  =  f0 
g  (qy)  H  g0 


d  t  = 


4> 


dx 

t  (x)  ^  9 

y  (x)  =  Cl  (I 


-  u)  F„  +  y 


(14) 


(15) 


(16) 


(17) 


Notice  that  Eq.  (17)  gives  the  exact  solution  of  the  linear  equation 
with  constant  coefficients. 

The  fundamental  question  to  be  explored  at  this  point,  the  trans¬ 
formation  having  been  obtained,  is  whether,  and  to  what  extent,  Eq.  (16) 
is  stiff.  Basic  to  this  question  is  the  matrix 


(- 


d  i, 
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which  is  associated  with  Eq,  (16)  in  the  same  manner  that  f  is  asso¬ 
ciated  with  Eq.  (2).  On  the  assumption  that  all  the  necessary  partial 
derivatives  exist,  it  can  be  shown  after  some  manipulation  that 


As 


m 


M 


it  follows  that  the  eigenvalues  of  m  are  the  same  as  the  eigenvalues 
of 


f  ( x ,  y )  -  f 0 


Thus,  the  spectral  radius  of  m  can  be  made  qs  small  as  desired 
simply  by  choosing  an  x  sufficiently  close  to  x0 .  It  is  assumed  here,  of 
course,  that  T  (x,  y)  is  continuous  at  x  =  x0 .  According  to  the  definition 
of  a  stiff  equation  given  earlier,  it  follows  that  Eq.  (16)  is  not  a  stiff 
equation  at  x  =  x0  +  h  provided  the  value  of  h  is  not  too  large.  Any 
failure  in  the  Runge-Kutta  integration  of  Eq.  (16)  will,  therefore,  be 
the  result  of  some  other  undesirable  phenomenon. 


SECTION  IV 

APPLICATION  TECHNIQUES 


One  method  of  applying  the  transformation  to  obtain  the  solution  of 
Eq.  (2)  would  be  to  solve  the  transformed  Eq.  (16)  for  t  (x)  by  the  Runge- 
Kutta  method,  compute  the  corresponding  y  (x) ,  update  the  transforma¬ 
tion,  and  continue  in  the  same  fashion  on  a  new  interval'.  While  t  (x) 
is  piecewise  discontinuous,  this  of  course  has  no  bearing  on  the  solution 
of  interest,  y  (x)  . 

Having  served  its  purpose,  the  transformation  can  now,  in  a  com¬ 
putational  sense,  be  entirely  removed  from  the  procedure.  Thus,  the 
technique  can  be  thought  of  in  a  completely  different  way,  that  of  solving 
the  untransformed  Eq.  (2)  by  a  method  different  from  that  of  Runge-Kutta. 

While  the  above  two  interpretations  of  the  method  would  theoretically 
yield  identical  answers,  the  second  version  seems  to  be  advantageous 
to  the  first  in  practice.  This  was  experienced  in  actual  test  runs  in 
which  the  two  interpretations  were  compared  on  the  basis  of  run  time, 
accuracy,  and  simplicity  of  program  logic. 
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As  a  result  of  Eqs.  (14)  and  (16),  and  the  classical  fourth-order 
Runge-Kutta  equations  (Ref,  6,  p.  122),  the  equations  necessary  for 
integrating  Eq.  (2)  with  the  transformation  removed  are 

for  j  -  0 

X  =  xc 


for  j  =  1 


for  j  =  2 


0 


x  =  x» +  ir 

Y  =  4"1  [i  -  5  toj  F„  +  yo 
kt  =  u“'  (x)^F  (x ,  y )  -  Fc  +  f0  (y  -  yo)J 

x  55  x°  +  T 

y  =  -i-  J  (x)  ^  +  f0~ 1  [V  -  u  (x)J  F0  +  yo 


K  =  5-‘  Wf?  (x,y)  -  F0  +  f„  (y  -  y#)l 

for  j  =  3 

x  =  x0  +  h 

y  =  h  U  (x)  kj  +  f~l  ["f  -  u  (X.]  F„  +  y0 

k3  =  u  “ '  (x)  £  F  ( x ,  y )  -  F0  +  T0  ( v  -  yo 
at 

x  =  x„  +  h 

y  =  ~  »  (*)  (  2k,  +  2ka  +  k} )  +  f [f  -  5  (x)J  Fq  +  y0 


(18a) 


(18b) 


(18c) 


(18  d) 


Notice  that  although  the  algebra  is  somewhat  more  involved,  there 
is  a  close  similarity  between\Eq.  (18)  and  the  corresponding  Runge- 
Kutta  equations. 


The  problem  that  initially  motivated  the  development  of  this  tech¬ 
nique  is  a  system  of  six  chemical  rate  equations  (Ref.  7).  However,  the 
matrix  T  (x,  y)  associated  with  this  problem  was  found  to  be  ill- 
conditioned;  hence,  a  modified  version  of  the  method  was  developed. 


Equation  (18)  in  scalar  form  was  applied  to  each  of  the  six  equa¬ 
tions  individually.  In  other  words,  the  six  equations  were,  for  purposes 
of  making  the  transformation  only,  assumed  to  be  uncoupled. 
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In  order  to  determine  whether  or  not  to  make  the  transformation  on 

the  i  th  equation,  — - f,(x,  y) ,  the  number  — ^  0  was  examined  to 

see  if  the  equation  was  suitable  for  integration  by  the  Runge-Kutta  Method. 

If 

-  -P-  Ax  >  0.5 
d  y 

then  the  i  th  equation  was  transformed.  If  the  above  inequality  did  not 
hold,  the  i  th  equation  was  considered  suitable  for  integration  without 
making  the  transformation  (Ref.  8,  p.  198,  or  Ref.  1).  Thus  the  six 
equations  were  solved  as  a  coupled  system,  but  the  integration  method 
used  on  the  individual  equations  was  dictated  by  the  above  criterion. 

Also,  the  test  was  made  at  every  integration  step  so  that  at  one  value 
of  x,  for  example,  Eqs.  (1),  (2),  (5),  and  (6)  were  transformed  while 
Eqs.  (2)  and  (4)  were  integrated  by  the  Runge-Kutta  method.  At  another 
value  of  x,  an  entirely  different  combination  might  very  well  have  pre¬ 
vailed. 


SECTION  V 

DISCUSSION  OF  RESULTS 


Several  stiff  systems  for  which  the  analytic  solution  is  known  were 
solved  both  by  the  Runge-Kutta  method  and  by  the  method  under  considera¬ 
tion  here.  One  such  equation  is 

y  '  +  a  y  =  xS  ;  y  (0)  =  0 

A  comparison  of  the  errors,  the  absolute  value  of  the  difference 
between  the  exact  answer  and  the  approximation,  is  shown  in  Fig.  1  for 
a  =  100  and  a  =  1000  . 

Figure  2  shows  a  similar  graph  for  the  Euclidian  norm  of  the  error 
vector  associated  with  the  system. 

-  y  '  +  a,,  yi  +  a12  y2  =  b,  x  ;  y;  (0)  =  jq® 

”  y'  +  a2.  >',  +  a22  y2  =  b2  xJ  ;  y2  (0)  =  y° 

No  analytic  solution  was  available  for  analyzing  the  success  of  the 
method  as  applied  to  the  chemical  rate  equations.  Essentially,  the  same 
solution  was  found  using  both  the  Runge-Kutta  method  proper  and  the 
transformed  method,  a  considerably  larger  Ax  being  admissible  in  the 
latter  case.  These  answers  were  also  compared  with  solutions  from 
another  source  (Ref.  7)  as  further  confirmation  of  their  validity. 
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SECTION  VI 
CONCLUSIONS 


The  method  discussed  herein  for  integrating  a  system  of  stiff,  first- 
order  differential  equations  has  been  found  to  be  practical  and  expedient 
in  all  cases  on  which  the  method  has  been  tested.  Answers  comparable 
to  those  indicated  by  the  Runge-Kutta  method  were  obtained  with  a  con¬ 
siderably  larger  integration  step  size. 
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